from sympy.integrals.laplace import (
    laplace_transform, inverse_laplace_transform,
    LaplaceTransform, InverseLaplaceTransform,
    _laplace_deep_collect, laplace_correspondence,
    laplace_initial_conds)
from sympy.core.function import Function, expand_mul
from sympy.core import EulerGamma, Subs, Derivative, diff
from sympy.core.exprtools import factor_terms
from sympy.core.numbers import I, oo, pi
from sympy.core.relational import Eq
from sympy.core.singleton import S
from sympy.core.symbol import Symbol, symbols
from sympy.simplify.simplify import simplify
from sympy.functions.elementary.complexes import Abs, re
from sympy.functions.elementary.exponential import exp, log, exp_polar
from sympy.functions.elementary.hyperbolic import cosh, sinh, coth, asinh
from sympy.functions.elementary.miscellaneous import sqrt
from sympy.functions.elementary.piecewise import Piecewise
from sympy.functions.elementary.trigonometric import atan, cos, sin
from sympy.functions.elementary.complexes import sign
from sympy.logic.boolalg import And
from sympy.functions.special.gamma_functions import (
    lowergamma, gamma, uppergamma)
from sympy.functions.special.delta_functions import DiracDelta, Heaviside
from sympy.functions.special.singularity_functions import SingularityFunction
from sympy.functions.special.zeta_functions import lerchphi
from sympy.functions.special.error_functions import (
    fresnelc, fresnels, erf, erfc, erfi, Ei, Ci, expint, E1)
from sympy.functions.special.bessel import besseli, besselj, besselk, bessely
from sympy.testing.pytest import slow, warns_deprecated_sympy, XFAIL
from sympy.matrices import Matrix, eye
from sympy.abc import s


@slow
def test_laplace_transform():
    LT = laplace_transform
    ILT = inverse_laplace_transform
    a, b, c = symbols('a, b, c', positive=True)
    np = symbols('np', integer=True, positive=True)
    t, w, x = symbols('t, w, x')
    f = Function('f')
    F = Function('F')
    g = Function('g')
    y = Function('y')
    Y = Function('Y')

    # Test helper functions
    assert (
        _laplace_deep_collect(exp((t+a)*(t+b)) +
                              besselj(2, exp((t+a)*(t+b)-t**2)), t) ==
        exp(a*b + t**2 + t*(a + b)) + besselj(2, exp(a*b + t*(a + b))))
    L = laplace_transform(diff(y(t), t, 3), t, s, noconds=True)
    L = laplace_correspondence(L, {y: Y})
    L = laplace_initial_conds(L, t, {y: [2, 4, 8, 16, 32]})
    assert L == s**3*Y(s) - 2*s**2 - 4*s - 8
    # Test whether `noconds=True` in `doit`:
    assert (2*LaplaceTransform(exp(t), t, s) - 1).doit() == -1 + 2/(s - 1)
    assert (LT(a*t+t**2+t**(S(5)/2), t, s) ==
            (a/s**2 + 2/s**3 + 15*sqrt(pi)/(8*s**(S(7)/2)), 0, True))
    assert LT(b/(t+a), t, s) == (-b*exp(-a*s)*Ei(-a*s), 0, True)
    assert (LT(1/sqrt(t+a), t, s) ==
            (sqrt(pi)*sqrt(1/s)*exp(a*s)*erfc(sqrt(a)*sqrt(s)), 0, True))
    assert (LT(sqrt(t)/(t+a), t, s) ==
            (-pi*sqrt(a)*exp(a*s)*erfc(sqrt(a)*sqrt(s)) + sqrt(pi)*sqrt(1/s),
             0, True))
    assert (LT((t+a)**(-S(3)/2), t, s) ==
            (-2*sqrt(pi)*sqrt(s)*exp(a*s)*erfc(sqrt(a)*sqrt(s)) + 2/sqrt(a),
             0, True))
    assert (LT(t**(S(1)/2)*(t+a)**(-1), t, s) ==
            (-pi*sqrt(a)*exp(a*s)*erfc(sqrt(a)*sqrt(s)) + sqrt(pi)*sqrt(1/s),
             0, True))
    assert (LT(1/(a*sqrt(t) + t**(3/2)), t, s) ==
            (pi*sqrt(a)*exp(a*s)*erfc(sqrt(a)*sqrt(s)), 0, True))
    assert (
        (LT((t+a)**b, t, s) ==
            (s**(-b - 1)*exp(a*s)*uppergamma(b + 1, a*s), 0, True))
        )
    assert LT(t**5/(t+a), t, s) == (120*a**5*uppergamma(-5, a*s), 0, True)
    assert LT(exp(t), t, s) == (1/(s - 1), 1, True)
    assert LT(exp(2*t), t, s) == (1/(s - 2), 2, True)
    assert LT(exp(a*t), t, s) == (1/(s - a), a, True)
    assert LT(exp(a*(t-b)), t, s) == (exp(-a*b)/(-a + s), a, True)
    assert LT(t*exp(-a*(t)), t, s) == ((a + s)**(-2), -a, True)
    assert LT(t*exp(-a*(t-b)), t, s) == (exp(a*b)/(a + s)**2, -a, True)
    assert LT(b*t*exp(-a*t), t, s) == (b/(a + s)**2, -a, True)
    assert LT(exp(-a*exp(-t)), t, s) == (lowergamma(s, a)/a**s, 0, True)
    assert LT(exp(-a*exp(t)), t, s) == (a**s*uppergamma(-s, a), 0, True)
    assert (LT(t**(S(7)/4)*exp(-8*t)/gamma(S(11)/4), t, s) ==
            ((s + 8)**(-S(11)/4), -8, True))
    assert (LT(t**(S(3)/2)*exp(-8*t), t, s) ==
            (3*sqrt(pi)/(4*(s + 8)**(S(5)/2)), -8, True))
    assert LT(t**a*exp(-a*t), t, s) == ((a+s)**(-a-1)*gamma(a+1), -a, True)
    assert (LT(b*exp(-a*t**2), t, s) ==
            (sqrt(pi)*b*exp(s**2/(4*a))*erfc(s/(2*sqrt(a)))/(2*sqrt(a)),
             0, True))
    assert (LT(exp(-2*t**2), t, s) ==
            (sqrt(2)*sqrt(pi)*exp(s**2/8)*erfc(sqrt(2)*s/4)/4, 0, True))
    # The following test case fails because when using the integral definition
    # of LT, a substitution into a divergent integral occurs, that does not
    # simplify into infinity. To be fixed in the future, for now it is marked
    # as an xfail at  test_laplace_transform_exp_2t2
    # assert (LT(b*exp(2*t**2), t, s) ==
    #         (b*LaplaceTransform(exp(2*t**2), t, s), -oo, True))
    assert (LT(t*exp(-a*t**2), t, s) ==
            (1/(2*a) - s*erfc(s/(2*sqrt(a)))/(4*sqrt(pi)*a**(S(3)/2)),
             0, True))
    assert (LT(exp(-a/t), t, s) ==
            (2*sqrt(a)*sqrt(1/s)*besselk(1, 2*sqrt(a)*sqrt(s)), 0, True))
    assert LT(sqrt(t)*exp(-a/t), t, s, simplify=True) == (
        sqrt(pi)*(sqrt(a)*sqrt(s) + 1/S(2))*sqrt(s**(-3)) *
        exp(-2*sqrt(a)*sqrt(s)), 0, True)
    assert (LT(exp(-a/t)/sqrt(t), t, s) ==
            (sqrt(pi)*sqrt(1/s)*exp(-2*sqrt(a)*sqrt(s)), 0, True))
    assert (LT(exp(-a/t)/(t*sqrt(t)), t, s) ==
            (sqrt(pi)*sqrt(1/a)*exp(-2*sqrt(a)*sqrt(s)), 0, True))
    # TODO: rules with sqrt(a*t) and sqrt(a/t) have stopped working after
    #       changes to as_base_exp
    # assert (
    #     LT(exp(-2*sqrt(a*t)), t, s) ==
    #     (1/s - sqrt(pi)*sqrt(a) * exp(a/s)*erfc(sqrt(a)*sqrt(1/s)) /
    #      s**(S(3)/2), 0, True))
    # assert LT(exp(-2*sqrt(a*t))/sqrt(t), t, s) == (
    #     exp(a/s)*erfc(sqrt(a) * sqrt(1/s))*(sqrt(pi)*sqrt(1/s)), 0, True)
    assert (LT(t**4*exp(-2/t), t, s) ==
            (8*sqrt(2)*(1/s)**(S(5)/2)*besselk(5, 2*sqrt(2)*sqrt(s)),
             0, True))
    assert LT(sinh(a*t), t, s) == (a/(-a**2 + s**2), a, True)
    assert (LT(b*sinh(a*t)**2, t, s) ==
            (2*a**2*b/(-4*a**2*s + s**3), 2*a, True))
    assert (LT(b*sinh(a*t)**2, t, s, simplify=True) ==
            (2*a**2*b/(s*(-4*a**2 + s**2)), 2*a, True))
    # The following line confirms that issue #21202 is solved
    assert LT(cosh(2*t), t, s) == (s/(-4 + s**2), 2, True)
    assert LT(cosh(a*t), t, s) == (s/(-a**2 + s**2), a, True)
    assert (LT(cosh(a*t)**2, t, s, simplify=True) ==
            ((2*a**2 - s**2)/(s*(4*a**2 - s**2)), 2*a, True))
    assert (LT(sinh(x+3), x, s, simplify=True) ==
            ((s*sinh(3) + cosh(3))/(s**2 - 1), 1, True))
    L, _, _ = LT(42*sin(w*t+x)**2, t, s)
    assert (
        L -
        21*(s**2 + s*(-s*cos(2*x) + 2*w*sin(2*x)) +
            4*w**2)/(s*(s**2 + 4*w**2))).simplify() == 0
    # The following line replaces the old test test_issue_7173()
    assert LT(sinh(a*t)*cosh(a*t), t, s, simplify=True) == (a/(-4*a**2 + s**2),
                                                            2*a, True)
    assert LT(sinh(a*t)/t, t, s) == (log((a + s)/(-a + s))/2, a, True)
    assert (LT(t**(-S(3)/2)*sinh(a*t), t, s) ==
            (-sqrt(pi)*(sqrt(-a + s) - sqrt(a + s)), a, True))
    # assert (LT(sinh(2*sqrt(a*t)), t, s) ==
    #         (sqrt(pi)*sqrt(a)*exp(a/s)/s**(S(3)/2), 0, True))
    # assert (LT(sqrt(t)*sinh(2*sqrt(a*t)), t, s, simplify=True) ==
    #         ((-sqrt(a)*s**(S(5)/2) + sqrt(pi)*s**2*(2*a + s)*exp(a/s) *
    #           erf(sqrt(a)*sqrt(1/s))/2)/s**(S(9)/2), 0, True))
    # assert (LT(sinh(2*sqrt(a*t))/sqrt(t), t, s) ==
    #         (sqrt(pi)*exp(a/s)*erf(sqrt(a)*sqrt(1/s))/sqrt(s), 0, True))
    # assert (LT(sinh(sqrt(a*t))**2/sqrt(t), t, s) ==
    #         (sqrt(pi)*(exp(a/s) - 1)/(2*sqrt(s)), 0, True))
    assert (LT(t**(S(3)/7)*cosh(a*t), t, s) ==
            (((a + s)**(-S(10)/7) + (-a+s)**(-S(10)/7))*gamma(S(10)/7)/2,
             a, True))
    # assert (LT(cosh(2*sqrt(a*t)), t, s) ==
    #         (sqrt(pi)*sqrt(a)*exp(a/s)*erf(sqrt(a)*sqrt(1/s))/s**(S(3)/2) +
    #          1/s, 0, True))
    # assert (LT(sqrt(t)*cosh(2*sqrt(a*t)), t, s) ==
    #         (sqrt(pi)*(a + s/2)*exp(a/s)/s**(S(5)/2), 0, True))
    # assert (LT(cosh(2*sqrt(a*t))/sqrt(t), t, s) ==
    #         (sqrt(pi)*exp(a/s)/sqrt(s), 0, True))
    # assert (LT(cosh(sqrt(a*t))**2/sqrt(t), t, s) ==
    #         (sqrt(pi)*(exp(a/s) + 1)/(2*sqrt(s)), 0, True))
    assert LT(log(t), t, s, simplify=True) == (
        (-log(s) - EulerGamma)/s, 0, True)
    assert (LT(-log(t/a), t, s, simplify=True) ==
            ((log(a) + log(s) + EulerGamma)/s, 0, True))
    assert LT(log(1+a*t), t, s) == (-exp(s/a)*Ei(-s/a)/s, 0, True)
    assert (LT(log(t+a), t, s, simplify=True) ==
            ((s*log(a) - exp(s/a)*Ei(-s/a))/s**2, 0, True))
    assert (LT(log(t)/sqrt(t), t, s, simplify=True) ==
            (sqrt(pi)*(-log(s) - log(4) - EulerGamma)/sqrt(s), 0, True))
    assert (LT(t**(S(5)/2)*log(t), t, s, simplify=True) ==
            (sqrt(pi)*(-15*log(s) - log(1073741824) - 15*EulerGamma + 46) /
             (8*s**(S(7)/2)), 0, True))
    assert (LT(t**3*log(t), t, s, noconds=True, simplify=True) -
            6*(-log(s) - S.EulerGamma + S(11)/6)/s**4).simplify() == S.Zero
    assert (LT(log(t)**2, t, s, simplify=True) ==
            (((log(s) + EulerGamma)**2 + pi**2/6)/s, 0, True))
    assert (LT(exp(-a*t)*log(t), t, s, simplify=True) ==
            ((-log(a + s) - EulerGamma)/(a + s), -a, True))
    assert LT(sin(a*t), t, s) == (a/(a**2 + s**2), 0, True)
    assert (LT(Abs(sin(a*t)), t, s) ==
            (a*coth(pi*s/(2*a))/(a**2 + s**2), 0, True))
    assert LT(sin(a*t)/t, t, s) == (atan(a/s), 0, True)
    assert LT(sin(a*t)**2/t, t, s) == (log(4*a**2/s**2 + 1)/4, 0, True)
    assert (LT(sin(a*t)**2/t**2, t, s) ==
            (a*atan(2*a/s) - s*log(4*a**2/s**2 + 1)/4, 0, True))
    # assert (LT(sin(2*sqrt(a*t)), t, s) ==
    #         (sqrt(pi)*sqrt(a)*exp(-a/s)/s**(S(3)/2), 0, True))
    # assert LT(sin(2*sqrt(a*t))/t, t, s) == (pi*erf(sqrt(a)*sqrt(1/s)), 0, True)
    assert LT(cos(a*t), t, s) == (s/(a**2 + s**2), 0, True)
    assert (LT(cos(a*t)**2, t, s) ==
            ((2*a**2 + s**2)/(s*(4*a**2 + s**2)), 0, True))
    # assert (LT(sqrt(t)*cos(2*sqrt(a*t)), t, s, simplify=True) ==
    #         (sqrt(pi)*(-a + s/2)*exp(-a/s)/s**(S(5)/2), 0, True))
    # assert (LT(cos(2*sqrt(a*t))/sqrt(t), t, s) ==
    #         (sqrt(pi)*sqrt(1/s)*exp(-a/s), 0, True))
    assert (LT(sin(a*t)*sin(b*t), t, s) ==
            (2*a*b*s/((s**2 + (a - b)**2)*(s**2 + (a + b)**2)), 0, True))
    assert (LT(cos(a*t)*sin(b*t), t, s) ==
            (b*(-a**2 + b**2 + s**2)/((s**2 + (a - b)**2)*(s**2 + (a + b)**2)),
             0, True))
    assert (LT(cos(a*t)*cos(b*t), t, s) ==
            (s*(a**2 + b**2 + s**2)/((s**2 + (a - b)**2)*(s**2 + (a + b)**2)),
             0, True))
    assert (LT(-a*t*cos(a*t) + sin(a*t), t, s, simplify=True) ==
            (2*a**3/(a**4 + 2*a**2*s**2 + s**4), 0, True))
    assert LT(c*exp(-b*t)*sin(a*t), t, s) == (a *
                                              c/(a**2 + (b + s)**2), -b, True)
    assert LT(c*exp(-b*t)*cos(a*t), t, s) == (c*(b + s)/(a**2 + (b + s)**2),
                                              -b, True)
    L, plane, cond = LT(cos(x + 3), x, s, simplify=True)
    assert plane == 0
    assert L - (s*cos(3) - sin(3))/(s**2 + 1) == 0
    # Error functions (laplace7.pdf)
    assert LT(erf(a*t), t, s) == (exp(s**2/(4*a**2))*erfc(s/(2*a))/s, 0, True)
    # assert LT(erf(sqrt(a*t)), t, s) == (sqrt(a)/(s*sqrt(a + s)), 0, True)
    # assert (LT(exp(a*t)*erf(sqrt(a*t)), t, s, simplify=True) ==
    #         (-sqrt(a)/(sqrt(s)*(a - s)), a, True))
    # assert (LT(erf(sqrt(a/t)/2), t, s, simplify=True) ==
    #         (1/s - exp(-sqrt(a)*sqrt(s))/s, 0, True))
    # assert (LT(erfc(sqrt(a*t)), t, s, simplify=True) ==
    #         (-sqrt(a)/(s*sqrt(a + s)) + 1/s, -a, True))
    # assert (LT(exp(a*t)*erfc(sqrt(a*t)), t, s) ==
    #         (1/(sqrt(a)*sqrt(s) + s), 0, True))
    # assert LT(erfc(sqrt(a/t)/2), t, s) == (exp(-sqrt(a)*sqrt(s))/s, 0, True)
    # Bessel functions (laplace8.pdf)
    assert LT(besselj(0, a*t), t, s) == (1/sqrt(a**2 + s**2), 0, True)
    assert (LT(besselj(1, a*t), t, s, simplify=True) ==
            (a/(a**2 + s**2 + s*sqrt(a**2 + s**2)), 0, True))
    assert (LT(besselj(2, a*t), t, s, simplify=True) ==
            (a**2/(sqrt(a**2 + s**2)*(s + sqrt(a**2 + s**2))**2), 0, True))
    assert (LT(t*besselj(0, a*t), t, s) ==
            (s/(a**2 + s**2)**(S(3)/2), 0, True))
    assert (LT(t*besselj(1, a*t), t, s) ==
            (a/(a**2 + s**2)**(S(3)/2), 0, True))
    assert (LT(t**2*besselj(2, a*t), t, s) ==
            (3*a**2/(a**2 + s**2)**(S(5)/2), 0, True))
    # assert LT(besselj(0, 2*sqrt(a*t)), t, s) == (exp(-a/s)/s, 0, True)
    # assert (LT(t**(S(3)/2)*besselj(3, 2*sqrt(a*t)), t, s) ==
    #         (a**(S(3)/2)*exp(-a/s)/s**4, 0, True))
    assert (LT(besselj(0, a*sqrt(t**2+b*t)), t, s, simplify=True) ==
            (exp(b*(s - sqrt(a**2 + s**2)))/sqrt(a**2 + s**2), 0, True))
    assert LT(besseli(0, a*t), t, s) == (1/sqrt(-a**2 + s**2), a, True)
    assert (LT(besseli(1, a*t), t, s, simplify=True) ==
            (a/(-a**2 + s**2 + s*sqrt(-a**2 + s**2)), a, True))
    assert (LT(besseli(2, a*t), t, s, simplify=True) ==
            (a**2/(sqrt(-a**2 + s**2)*(s + sqrt(-a**2 + s**2))**2), a, True))
    assert LT(t*besseli(0, a*t), t, s) == (s/(-a**2 + s**2)**(S(3)/2), a, True)
    assert LT(t*besseli(1, a*t), t, s) == (a/(-a**2 + s**2)**(S(3)/2), a, True)
    assert (LT(t**2*besseli(2, a*t), t, s) ==
            (3*a**2/(-a**2 + s**2)**(S(5)/2), a, True))
    # assert (LT(t**(S(3)/2)*besseli(3, 2*sqrt(a*t)), t, s) ==
    #         (a**(S(3)/2)*exp(a/s)/s**4, 0, True))
    assert (LT(bessely(0, a*t), t, s) ==
            (-2*asinh(s/a)/(pi*sqrt(a**2 + s**2)), 0, True))
    assert (LT(besselk(0, a*t), t, s) ==
            (log((s + sqrt(-a**2 + s**2))/a)/sqrt(-a**2 + s**2), -a, True))
    assert (LT(sin(a*t)**4, t, s, simplify=True) ==
            (24*a**4/(s*(64*a**4 + 20*a**2*s**2 + s**4)), 0, True))
    # Test general rules and unevaluated forms
    # These all also test whether issue #7219 is solved.
    assert LT(Heaviside(t-1)*cos(t-1), t, s) == (s*exp(-s)/(s**2 + 1), 0, True)
    assert LT(a*f(t), t, w) == (a*LaplaceTransform(f(t), t, w), -oo, True)
    assert (LT(a*Heaviside(t+1)*f(t+1), t, s) ==
            (a*LaplaceTransform(f(t + 1), t, s), -oo, True))
    assert (LT(a*Heaviside(t-1)*f(t-1), t, s) ==
            (a*LaplaceTransform(f(t), t, s)*exp(-s), -oo, True))
    assert (LT(b*f(t/a), t, s) ==
            (a*b*LaplaceTransform(f(t), t, a*s), -oo, True))
    assert LT(exp(-f(x)*t), t, s) == (1/(s + f(x)), -re(f(x)), True)
    assert (LT(exp(-a*t)*f(t), t, s) ==
            (LaplaceTransform(f(t), t, a + s), -oo, True))
    # assert (LT(exp(-a*t)*erfc(sqrt(b/t)/2), t, s) ==
    #         (exp(-sqrt(b)*sqrt(a + s))/(a + s), -a, True))
    assert (LT(sinh(a*t)*f(t), t, s) ==
            (LaplaceTransform(f(t), t, -a + s)/2 -
             LaplaceTransform(f(t), t, a + s)/2, -oo, True))
    assert (LT(sinh(a*t)*t, t, s, simplify=True) ==
            (2*a*s/(a**4 - 2*a**2*s**2 + s**4), a, True))
    assert (LT(cosh(a*t)*f(t), t, s) ==
            (LaplaceTransform(f(t), t, -a + s)/2 +
             LaplaceTransform(f(t), t, a + s)/2, -oo, True))
    assert (LT(cosh(a*t)*t, t, s, simplify=True) ==
            (1/(2*(a + s)**2) + 1/(2*(a - s)**2), a, True))
    assert (LT(sin(a*t)*f(t), t, s, simplify=True) ==
            (I*(-LaplaceTransform(f(t), t, -I*a + s) +
                LaplaceTransform(f(t), t, I*a + s))/2, -oo, True))
    assert (LT(sin(f(t)), t, s) ==
            (LaplaceTransform(sin(f(t)), t, s), -oo, True))
    assert (LT(sin(a*t)*t, t, s, simplify=True) ==
            (2*a*s/(a**4 + 2*a**2*s**2 + s**4), 0, True))
    assert (LT(cos(a*t)*f(t), t, s) ==
            (LaplaceTransform(f(t), t, -I*a + s)/2 +
             LaplaceTransform(f(t), t, I*a + s)/2, -oo, True))
    assert (LT(cos(a*t)*t, t, s, simplify=True) ==
            ((-a**2 + s**2)/(a**4 + 2*a**2*s**2 + s**4), 0, True))
    L, plane, _ = LT(sin(a*t+b)**2*f(t), t, s)
    assert plane == -oo
    assert (
        -L + (
            LaplaceTransform(f(t), t, s)/2 -
            LaplaceTransform(f(t), t, -2*I*a + s)*exp(2*I*b)/4 -
            LaplaceTransform(f(t), t, 2*I*a + s)*exp(-2*I*b)/4)) == 0
    L = LT(sin(a*t+b)**2*f(t), t, s, noconds=True)
    assert (
        laplace_correspondence(L, {f: F}) ==
        F(s)/2 - F(-2*I*a + s)*exp(2*I*b)/4 -
        F(2*I*a + s)*exp(-2*I*b)/4)
    L, plane, _ = LT(sin(a*t)**3*cosh(b*t), t, s)
    assert plane == b
    assert (
        -L - 3*a/(8*(9*a**2 + b**2 + 2*b*s + s**2)) -
        3*a/(8*(9*a**2 + b**2 - 2*b*s + s**2)) +
        3*a/(8*(a**2 + b**2 + 2*b*s + s**2)) +
        3*a/(8*(a**2 + b**2 - 2*b*s + s**2))).simplify() == 0
    assert (LT(t**2*exp(-t**2), t, s) ==
            (sqrt(pi)*s**2*exp(s**2/4)*erfc(s/2)/8 - s/4 +
             sqrt(pi)*exp(s**2/4)*erfc(s/2)/4, 0, True))
    assert (LT((a*t**2 + b*t + c)*f(t), t, s) ==
            (a*Derivative(LaplaceTransform(f(t), t, s), (s, 2)) -
             b*Derivative(LaplaceTransform(f(t), t, s), s) +
            c*LaplaceTransform(f(t), t, s), -oo, True))
    assert (LT(t**np*g(t), t, s) ==
            ((-1)**np*Derivative(LaplaceTransform(g(t), t, s), (s, np)),
             -oo, True))
    # The following tests check whether _piecewise_to_heaviside works:
    x1 = Piecewise((0, t <= 0), (1, t <= 1), (0, True))
    X1 = LT(x1, t, s)[0]
    assert X1 == 1/s - exp(-s)/s
    y1 = ILT(X1, s, t)
    assert y1 == Heaviside(t) - Heaviside(t - 1)
    x1 = Piecewise((0, t <= 0), (t, t <= 1), (2-t, t <= 2), (0, True))
    X1 = LT(x1, t, s)[0].simplify()
    assert X1 == (exp(2*s) - 2*exp(s) + 1)*exp(-2*s)/s**2
    y1 = ILT(X1, s, t)
    assert (
        -y1 + t*Heaviside(t) + (t - 2)*Heaviside(t - 2) -
        2*(t - 1)*Heaviside(t - 1)).simplify() == 0
    x1 = Piecewise((exp(t), t <= 0), (1, t <= 1), (exp(-(t)), True))
    X1 = LT(x1, t, s)[0]
    assert X1 == exp(-1)*exp(-s)/(s + 1) + 1/s - exp(-s)/s
    y1 = ILT(X1, s, t)
    assert y1 == (
        exp(-1)*exp(1 - t)*Heaviside(t - 1) + Heaviside(t) - Heaviside(t - 1))
    x1 = Piecewise((0, x <= 0), (1, x <= 1), (0, True))
    X1 = LT(x1, t, s)[0]
    assert X1 == Piecewise((0, x <= 0), (1, x <= 1), (0, True))/s
    x1 = [
        a*Piecewise((1, And(t > 1, t <= 3)), (2, True)),
        a*Piecewise((1, And(t >= 1, t <= 3)), (2, True)),
        a*Piecewise((1, And(t >= 1, t < 3)), (2, True)),
        a*Piecewise((1, And(t > 1, t < 3)), (2, True))]
    for x2 in x1:
        assert LT(x2, t, s)[0].expand() == 2*a/s - a*exp(-s)/s + a*exp(-3*s)/s
    assert (
        LT(Piecewise((1, Eq(t, 1)), (2, True)), t, s)[0] ==
        LaplaceTransform(Piecewise((1, Eq(t, 1)), (2, True)), t, s))
    # The following lines test whether _laplace_transform successfully
    # removes Heaviside(1) before processing espressions. It fails if
    # Heaviside(t) remains because then meijerg functions will appear.
    X1 = 1/sqrt(a*s**2-b)
    x1 = ILT(X1, s, t)
    Y1 = LT(x1, t, s)[0]
    Z1 = (Y1**2/X1**2).simplify()
    assert Z1 == 1
    # The following two lines test whether issues #5813 and #7176 are solved.
    assert (LT(diff(f(t), (t, 1)), t, s, noconds=True) ==
            s*LaplaceTransform(f(t), t, s) - f(0))
    assert (LT(diff(f(t), (t, 3)), t, s, noconds=True) ==
            s**3*LaplaceTransform(f(t), t, s) - s**2*f(0) -
            s*Subs(Derivative(f(t), t), t, 0) -
            Subs(Derivative(f(t), (t, 2)), t, 0))
    # Issue #7219
    assert (LT(diff(f(x, t, w), t, 2), t, s) ==
            (s**2*LaplaceTransform(f(x, t, w), t, s) - s*f(x, 0, w) -
             Subs(Derivative(f(x, t, w), t), t, 0), -oo, True))
    # Issue #23307
    assert (LT(10*diff(f(t), (t, 1)), t, s, noconds=True) ==
            10*s*LaplaceTransform(f(t), t, s) - 10*f(0))
    assert (LT(a*f(b*t)+g(c*t), t, s, noconds=True) ==
            a*LaplaceTransform(f(t), t, s/b)/b +
            LaplaceTransform(g(t), t, s/c)/c)
    assert inverse_laplace_transform(
        f(w), w, t, plane=0) == InverseLaplaceTransform(f(w), w, t, 0)
    assert (LT(f(t)*g(t), t, s, noconds=True) ==
            LaplaceTransform(f(t)*g(t), t, s))
    # Issue #24294
    assert (LT(b*f(a*t), t, s, noconds=True) ==
            b*LaplaceTransform(f(t), t, s/a)/a)
    assert LT(3*exp(t)*Heaviside(t), t, s) == (3/(s - 1), 1, True)
    assert (LT(2*sin(t)*Heaviside(t), t, s, simplify=True) ==
            (2/(s**2 + 1), 0, True))
    # Issue #25293
    assert (
        LT((1/(t-1))*sin(4*pi*(t-1))*DiracDelta(t-1) *
           (Heaviside(t-1/4) - Heaviside(t-2)), t, s)[0] == 4*pi*exp(-s))
    # Issue #28686
    assert (LT((t + 4)*a, t, s) == (a*(4*s + 1)/s**2, 0, True))
    # additional basic tests from wikipedia
    assert (LT((t - a)**b*exp(-c*(t - a))*Heaviside(t - a), t, s) ==
            ((c + s)**(-b - 1)*exp(-a*s)*gamma(b + 1), -c, True))
    assert (
        LT((exp(2*t)-1)*exp(-b-t)*Heaviside(t)/2, t, s, noconds=True,
           simplify=True) ==
        exp(-b)/(s**2 - 1))
    # DiracDelta function: standard cases
    assert LT(DiracDelta(t), t, s) == (1, -oo, True)
    assert LT(DiracDelta(a*t), t, s) == (1/a, -oo, True)
    assert LT(DiracDelta(t/42), t, s) == (42, -oo, True)
    assert LT(DiracDelta(t+42), t, s) == (0, -oo, True)
    assert (LT(DiracDelta(t)+DiracDelta(t-42), t, s) ==
            (1 + exp(-42*s), -oo, True))
    assert (LT(DiracDelta(t)-a*exp(-a*t), t, s, simplify=True) ==
            (s/(a + s), -a, True))
    assert (
        LT(exp(-t)*(DiracDelta(t)+DiracDelta(t-42)), t, s, simplify=True) ==
        (exp(-42*s - 42) + 1, -oo, True))
    assert LT(f(t)*DiracDelta(t-42), t, s) == (f(42)*exp(-42*s), -oo, True)
    assert LT(f(t)*DiracDelta(b*t-a), t, s) == (f(a/b)*exp(-a*s/b)/b,
                                                -oo, True)
    assert LT(f(t)*DiracDelta(b*t+a), t, s) == (0, -oo, True)
    # SingularityFunction
    assert LT(SingularityFunction(t, a, -1), t, s)[0] == exp(-a*s)
    assert LT(SingularityFunction(t, a, 1), t, s)[0] == exp(-a*s)/s**2
    assert LT(SingularityFunction(t, a, x), t, s)[0] == (
        LaplaceTransform(SingularityFunction(t, a, x), t, s))
    # Collection of cases that cannot be fully evaluated and/or would catch
    # some common implementation errors
    assert (LT(DiracDelta(t**2), t, s, noconds=True) ==
            LaplaceTransform(DiracDelta(t**2), t, s))
    assert LT(DiracDelta(t**2 - 1), t, s) == (exp(-s)/2, -oo, True)
    assert LT(DiracDelta(t*(1 - t)), t, s) == (1 - exp(-s), -oo, True)
    assert (LT((DiracDelta(t) + 1)*(DiracDelta(t - 1) + 1), t, s) ==
            (LaplaceTransform(DiracDelta(t)*DiracDelta(t - 1), t, s) +
             1 + exp(-s) + 1/s, 0, True))
    assert LT(DiracDelta(2*t-2*exp(a)), t, s) == (exp(-s*exp(a))/2, -oo, True)
    assert LT(DiracDelta(-2*t+2*exp(a)), t, s) == (exp(-s*exp(a))/2, -oo, True)
    # Heaviside tests
    assert LT(Heaviside(t), t, s) == (1/s, 0, True)
    assert LT(Heaviside(t - a), t, s) == (exp(-a*s)/s, 0, True)
    assert LT(Heaviside(t-1), t, s) == (exp(-s)/s, 0, True)
    assert LT(Heaviside(2*t-4), t, s) == (exp(-2*s)/s, 0, True)
    assert LT(Heaviside(2*t+4), t, s) == (1/s, 0, True)
    assert (LT(Heaviside(-2*t+4), t, s, simplify=True) ==
            (1/s - exp(-2*s)/s, 0, True))
    assert (LT(g(t)*Heaviside(t - w), t, s) ==
            (LaplaceTransform(g(t)*Heaviside(t - w), t, s), -oo, True))
    assert (
        LT(Heaviside(t-a)*g(t), t, s) ==
        (LaplaceTransform(g(a + t), t, s)*exp(-a*s), -oo, True))
    assert (
        LT(Heaviside(t+a)*g(t), t, s) ==
        (LaplaceTransform(g(t), t, s), -oo, True))
    assert (
        LT(Heaviside(-t+a)*g(t), t, s) ==
        (LaplaceTransform(g(t), t, s) -
         LaplaceTransform(g(a + t), t, s)*exp(-a*s), -oo, True))
    assert (
        LT(Heaviside(-t-a)*g(t), t, s) == (0, 0, True))
    # Fresnel functions
    assert (laplace_transform(fresnels(t), t, s, simplify=True) ==
            ((-sin(s**2/(2*pi))*fresnels(s/pi) +
              sqrt(2)*sin(s**2/(2*pi) + pi/4)/2 -
              cos(s**2/(2*pi))*fresnelc(s/pi))/s, 0, True))
    assert (laplace_transform(fresnelc(t), t, s, simplify=True) ==
            ((sin(s**2/(2*pi))*fresnelc(s/pi) -
              cos(s**2/(2*pi))*fresnels(s/pi) +
              sqrt(2)*cos(s**2/(2*pi) + pi/4)/2)/s, 0, True))
    # Matrix tests
    Mt = Matrix([[exp(t), t*exp(-t)], [t*exp(-t), exp(t)]])
    Ms = Matrix([[1/(s - 1), (s + 1)**(-2)],
                 [(s + 1)**(-2),     1/(s - 1)]])
    # The default behaviour for Laplace transform of a Matrix returns a Matrix
    # of Tuples and is deprecated:
    with warns_deprecated_sympy():
        Ms_conds = Matrix(
            [[(1/(s - 1), 1, True), ((s + 1)**(-2), -1, True)],
             [((s + 1)**(-2), -1, True), (1/(s - 1), 1, True)]])
    with warns_deprecated_sympy():
        assert LT(Mt, t, s) == Ms_conds
    # The new behavior is to return a tuple of a Matrix and the convergence
    # conditions for the matrix as a whole:
    assert LT(Mt, t, s, legacy_matrix=False) == (Ms, 1, True)
    # With noconds=True the transformed matrix is returned without conditions
    # either way:
    assert LT(Mt, t, s, noconds=True) == Ms
    assert LT(Mt, t, s, legacy_matrix=False, noconds=True) == Ms


@slow
@XFAIL
def test_laplace_transform_exp_2t2():
    # the test belongs in test_laplace_transform once fixed and is already
    # there but commented out
    LT = laplace_transform
    b = symbols('b', positive=True)
    t = symbols('t')
    wrong = (b*(oo*sign(exp(-s**2/8)) + sqrt(2)*sqrt(pi)*exp(-s**2/8)*erfi(sqrt(2)*s/4)/4), -oo, True)
    assert (LT(b*exp(2*t**2), t, s) != wrong)


@slow
def test_inverse_laplace_transform():
    s = symbols('s')
    k, n, t = symbols('k, n, t', real=True)
    a, b, c, d = symbols('a, b, c, d', positive=True)
    f = Function('f')
    F = Function('F')

    def ILT(g):
        return inverse_laplace_transform(g, s, t)

    def ILTS(g):
        return inverse_laplace_transform(g, s, t, simplify=True)

    def ILTF(g):
        return laplace_correspondence(
            inverse_laplace_transform(g, s, t), {f: F})

    # Tests for the rules in Bateman54.

    # Section 4.1: Some of the Laplace transform rules can also be used well
    #     in the inverse transform.
    assert ILTF(exp(-a*s)*F(s)) == f(-a + t)
    assert ILTF(k*F(s-a)) == k*f(t)*exp(-a*t)
    assert ILTF(diff(F(s), s, 3)) == -t**3*f(t)
    assert ILTF(diff(F(s), s, 4)) == t**4*f(t)

    # Section 5.1: Most rules are impractical for a computer algebra system.

    # Section 5.2: Rational functions
    assert ILT(2) == 2*DiracDelta(t)
    assert ILT(1/s) == Heaviside(t)
    assert ILT(1/s**2) == t*Heaviside(t)
    assert ILT(1/s**5) == t**4*Heaviside(t)/24
    assert ILT(1/s**n) == t**(n - 1)*Heaviside(t)/gamma(n)
    assert ILT(a/(a + s)) == a*exp(-a*t)*Heaviside(t)
    assert ILT(s/(a + s)) == -a*exp(-a*t)*Heaviside(t) + DiracDelta(t)
    assert (ILT(b*s/(s+a)**2) ==
            b*(-a*t*exp(-a*t)*Heaviside(t) + exp(-a*t)*Heaviside(t)))
    assert (ILTS(c/((s+a)*(s+b))) ==
            c*(exp(a*t) - exp(b*t))*exp(-t*(a + b))*Heaviside(t)/(a - b))
    assert (ILTS(c*s/((s+a)*(s+b))) ==
            c*(a*exp(b*t) - b*exp(a*t))*exp(-t*(a + b))*Heaviside(t)/(a - b))
    assert ILTS(s/(a + s)**3) == t*(-a*t + 2)*exp(-a*t)*Heaviside(t)/2
    assert ILTS(1/(s*(a + s)**3)) == (
        -a**2*t**2 - 2*a*t + 2*exp(a*t) - 2)*exp(-a*t)*Heaviside(t)/(2*a**3)
    assert ILT(1/(s*(a + s)**n)) == (
        Heaviside(t)*lowergamma(n, a*t)/(a**n*gamma(n)))
    assert ILT((s-a)**(-b)) == t**(b - 1)*exp(a*t)*Heaviside(t)/gamma(b)
    assert ILT((a + s)**(-2)) == t*exp(-a*t)*Heaviside(t)
    assert ILT((a + s)**(-5)) == t**4*exp(-a*t)*Heaviside(t)/24
    assert ILT(s**2/(s**2 + 1)) == -sin(t)*Heaviside(t) + DiracDelta(t)
    assert ILT(1 - 1/(s**2 + 1)) == -sin(t)*Heaviside(t) + DiracDelta(t)
    assert ILT(a/(a**2 + s**2)) == sin(a*t)*Heaviside(t)
    assert ILT(s/(s**2 + a**2)) == cos(a*t)*Heaviside(t)
    assert ILT(b/(b**2 + (a + s)**2)) == exp(-a*t)*sin(b*t)*Heaviside(t)
    assert (ILT(b*s/(b**2 + (a + s)**2)) ==
            b*(-a*exp(-a*t)*sin(b*t)/b + exp(-a*t)*cos(b*t))*Heaviside(t))
    assert ILT(1/(s**2*(s**2 + 1))) == t*Heaviside(t) - sin(t)*Heaviside(t)
    assert (ILTS(c*s/(d**2*(s+a)**2+b**2)) ==
            c*(-a*d*sin(b*t/d) + b*cos(b*t/d))*exp(-a*t)*Heaviside(t)/(b*d**2))
    assert ILTS((b*s**2 + d)/(a**2 + s**2)**2) == (
        2*a**2*b*sin(a*t) + (a**2*b - d)*(a*t*cos(a*t) -
                                          sin(a*t)))*Heaviside(t)/(2*a**3)
    assert ILTS(b/(s**2-a**2)) == b*sinh(a*t)*Heaviside(t)/a
    assert (ILT(b/(s**2-a**2)) ==
            b*(exp(a*t)*Heaviside(t)/(2*a) - exp(-a*t)*Heaviside(t)/(2*a)))
    assert ILTS(b*s/(s**2-a**2)) == b*cosh(a*t)*Heaviside(t)
    assert (ILT(b/(s*(s+a))) ==
            b*(Heaviside(t)/a - exp(-a*t)*Heaviside(t)/a))
    # Issue #24424
    assert (ILTS((s + 8)/((s + 2)*(s**2 + 2*s + 10))) ==
            ((8*sin(3*t) - 9*cos(3*t))*exp(t) + 9)*exp(-2*t)*Heaviside(t)/15)
    # Issue #8514; this is not important anymore, since this function
    # is not solved by integration anymore
    assert (ILT(1/(a*s**2+b*s+c)) ==
            2*exp(-b*t/(2*a))*sin(t*sqrt(4*a*c - b**2)/(2*a)) *
            Heaviside(t)/sqrt(4*a*c - b**2))

    # Section 5.3: Irrational algebraic functions
    assert (  # (1)
        ILT(1/sqrt(s)/(b*s-a)) ==
        exp(a*t/b)*Heaviside(t)*erf(sqrt(a)*sqrt(t)/sqrt(b))/(sqrt(a)*sqrt(b)))
    assert (  # (2)
        ILT(1/sqrt(k*s)/(c*s-a)/s) ==
        (-2*c*sqrt(t)/(sqrt(pi)*a) +
         c**(S(3)/2)*exp(a*t/c)*erf(sqrt(a)*sqrt(t)/sqrt(c))/a**(S(3)/2)) *
        Heaviside(t)/(c*sqrt(k)))
    assert (  # (4)
        ILT(1/(sqrt(c*s)+a)) == (-a*exp(a**2*t/c)*erfc(a*sqrt(t)/sqrt(c))/c +
                                 1/(sqrt(pi)*sqrt(c)*sqrt(t)))*Heaviside(t))
    assert (  # (5)
        ILT(a/s/(b*sqrt(s)+a)) ==
        (-exp(a**2*t/b**2)*erfc(a*sqrt(t)/b) + 1)*Heaviside(t))
    assert (  # (6)
            ILT((a-b)*sqrt(s)/(sqrt(s)+sqrt(a))/(s-b)) ==
            (sqrt(a)*sqrt(b)*exp(b*t)*erfc(sqrt(b)*sqrt(t)) +
             a*exp(a*t)*erfc(sqrt(a)*sqrt(t)) - b*exp(b*t))*Heaviside(t))
    assert (  # (7)
            ILT(1/sqrt(s)/(sqrt(b*s)+a)) ==
            exp(a**2*t/b)*Heaviside(t)*erfc(a*sqrt(t)/sqrt(b))/sqrt(b))
    assert (  # (8)
            ILT(a**2/(sqrt(s)+a)/s**(S(3)/2)) ==
            (2*a*sqrt(t)/sqrt(pi) + exp(a**2*t)*erfc(a*sqrt(t)) - 1) *
            Heaviside(t))
    assert (  # (9)
            ILT((a-b)*sqrt(b)/(s-b)/sqrt(s)/(sqrt(s)+sqrt(a))) ==
            (sqrt(a)*exp(b*t)*erf(sqrt(b)*sqrt(t)) +
             sqrt(b)*exp(a*t)*erfc(sqrt(a)*sqrt(t)) -
             sqrt(b)*exp(b*t))*Heaviside(t))
    assert (  # (10)
            ILT(1/(sqrt(s)+sqrt(a))**2) ==
            (-2*sqrt(a)*sqrt(t)/sqrt(pi) +
             (-2*a*t + 1)*(erf(sqrt(a)*sqrt(t)) -
                           1)*exp(a*t) + 1)*Heaviside(t))
    assert (  # (11)
            ILT(1/(sqrt(s)+sqrt(a))**2/s) ==
            ((2*t - 1/a)*exp(a*t)*erfc(sqrt(a)*sqrt(t)) + 1/a -
             2*sqrt(t)/(sqrt(pi)*sqrt(a)))*Heaviside(t))
    assert (  # (12)
            ILT(1/(sqrt(s)+a)**2/sqrt(s)) ==
            (-2*a*t*exp(a**2*t)*erfc(a*sqrt(t)) +
             2*sqrt(t)/sqrt(pi))*Heaviside(t))
    assert (  # (13)
            ILT(1/(sqrt(s)+a)**3) ==
            (-a*t*(2*a**2*t + 3)*exp(a**2*t)*erfc(a*sqrt(t)) +
             2*sqrt(t)*(a**2*t + 1)/sqrt(pi))*Heaviside(t))
    x = (
        - ILT(sqrt(s)/(sqrt(s)+a)**3) +
        2*(sqrt(pi)*a**2*t*(-2*sqrt(pi)*erfc(a*sqrt(t)) +
                            2*exp(-a**2*t)/(a*sqrt(t))) *
           (-a**4*t**2 - 5*a**2*t/2 - S.Half) * exp(a**2*t)/2 +
           sqrt(pi)*a*sqrt(t)*(a**2*t + 1)/2) *
        Heaviside(t)/(pi*a**2*t)).simplify()
    assert (  # (14)
            x == 0)
    x = (
        - ILT(1/sqrt(s)/(sqrt(s)+a)**3) +
        Heaviside(t)*(sqrt(t)*((2*a**2*t + 1) *
                               (sqrt(pi)*a*sqrt(t)*exp(a**2*t) *
                                erfc(a*sqrt(t)) - 1) + 1) /
                      (sqrt(pi)*a))).simplify()
    assert (  # (15)
            x == 0)
    assert (  # (16)
            factor_terms(ILT(3/(sqrt(s)+a)**4)) ==
            3*(-2*a**3*t**(S(5)/2)*(2*a**2*t + 5)/(3*sqrt(pi)) +
               t*(4*a**4*t**2 + 12*a**2*t + 3)*exp(a**2*t) *
               erfc(a*sqrt(t))/3)*Heaviside(t))
    assert (  # (17)
            ILT((sqrt(s)-a)/(s*(sqrt(s)+a))) ==
            (2*exp(a**2*t)*erfc(a*sqrt(t))-1)*Heaviside(t))
    assert (  # (18)
            ILT((sqrt(s)-a)**2/(s*(sqrt(s)+a)**2)) == (
                1 + 8*a**2*t*exp(a**2*t)*erfc(a*sqrt(t)) -
                8/sqrt(pi)*a*sqrt(t))*Heaviside(t))
    assert (  # (19)
            ILT((sqrt(s)-a)**3/(s*(sqrt(s)+a)**3)) == Heaviside(t)*(
                2*(8*a**4*t**2+8*a**2*t+1)*exp(a**2*t) *
                erfc(a*sqrt(t))-8/sqrt(pi)*a*sqrt(t)*(2*a**2*t+1)-1))
    assert (  # (22)
            ILT(sqrt(s+a)/(s+b)) == Heaviside(t)*(
                exp(-a*t)/sqrt(t)/sqrt(pi) +
                sqrt(a-b)*exp(-b*t)*erf(sqrt(a-b)*sqrt(t))))
    assert (  # (23)
            ILT(1/sqrt(s+b)/(s+a)) == Heaviside(t)*(
                1/sqrt(b-a)*exp(-a*t)*erf(sqrt(b-a)*sqrt(t))))
    assert (  # (35)
            ILT(1/sqrt(s**2+a**2)) == Heaviside(t)*(
                besselj(0, a*t)))
    assert (  # (44)
            ILT(1/sqrt(s**2-a**2)) == Heaviside(t)*(
                besseli(0, a*t)))

    # Miscellaneous tests
    # Can _inverse_laplace_time_shift deal with positive exponents?
    assert (
        - ILT((s**2*exp(2*s) + 4*exp(s) - 4)*exp(-2*s)/(s*(s**2 + 1))) +
        cos(t)*Heaviside(t) + 4*cos(t - 2)*Heaviside(t - 2) -
        4*cos(t - 1)*Heaviside(t - 1) - 4*Heaviside(t - 2) +
        4*Heaviside(t - 1)).simplify() == 0


@slow
def test_inverse_laplace_transform_old():
    from sympy.functions.special.delta_functions import DiracDelta
    ILT = inverse_laplace_transform
    a, b, c, d = symbols('a b c d', positive=True)
    n, r = symbols('n, r', real=True)
    t, z = symbols('t z')
    f = Function('f')
    F = Function('F')

    def simp_hyp(expr):
        return factor_terms(expand_mul(expr)).rewrite(sin)

    L = ILT(F(s), s, t)
    assert laplace_correspondence(L, {f: F}) == f(t)
    assert ILT(exp(-a*s)/s, s, t) == Heaviside(-a + t)
    assert ILT(exp(-a*s)/(b + s), s, t) == exp(-b*(-a + t))*Heaviside(-a + t)
    assert (ILT((b + s)/(a**2 + (b + s)**2), s, t) ==
            exp(-b*t)*cos(a*t)*Heaviside(t))
    assert (ILT(exp(-a*s)/s**b, s, t) ==
            (-a + t)**(b - 1)*Heaviside(-a + t)/gamma(b))
    assert (ILT(exp(-a*s)/sqrt(s**2 + 1), s, t) ==
            Heaviside(-a + t)*besselj(0, a - t))
    assert ILT(1/(s*sqrt(s + 1)), s, t) == Heaviside(t)*erf(sqrt(t))
    # TODO sinh/cosh shifted come out a mess. also delayed trig is a mess
    # TODO should this simplify further?
    assert (ILT(exp(-a*s)/s**b, s, t) ==
            (t - a)**(b - 1)*Heaviside(t - a)/gamma(b))
    assert (ILT(exp(-a*s)/sqrt(1 + s**2), s, t) ==
            Heaviside(t - a)*besselj(0, a - t))  # note: besselj(0, x) is even
    # XXX ILT turns these branch factor into trig functions ...
    assert (
        simplify(ILT(a**b*(s + sqrt(s**2 - a**2))**(-b)/sqrt(s**2 - a**2),
                     s, t).rewrite(exp)) ==
        Heaviside(t)*besseli(b, a*t))
    assert (
        ILT(a**b*(s + sqrt(s**2 + a**2))**(-b)/sqrt(s**2 + a**2),
            s, t, simplify=True).rewrite(exp) ==
        Heaviside(t)*besselj(b, a*t))
    assert ILT(1/(s*sqrt(s + 1)), s, t) == Heaviside(t)*erf(sqrt(t))
    # TODO can we make erf(t) work?
    assert (ILT((s * eye(2) - Matrix([[1, 0], [0, 2]])).inv(), s, t) ==
            Matrix([[exp(t)*Heaviside(t), 0], [0, exp(2*t)*Heaviside(t)]]))
    # Test time_diff rule
    assert (ILT(s**42*f(s), s, t) ==
            Derivative(InverseLaplaceTransform(f(s), s, t, None), (t, 42)))
    assert ILT(cos(s), s, t) == InverseLaplaceTransform(cos(s), s, t, None)
    # Rules for testing different DiracDelta cases
    assert (
        ILT(1 + 2*s + 3*s**2 + 5*s**3, s, t) == DiracDelta(t) +
        2*DiracDelta(t, 1) + 3*DiracDelta(t, 2) + 5*DiracDelta(t, 3))
    assert (ILT(2*exp(3*s) - 5*exp(-7*s), s, t) ==
            2*InverseLaplaceTransform(exp(3*s), s, t, None) -
            5*DiracDelta(t - 7))
    a = cos(sin(7)/2)
    assert ILT(a*exp(-3*s), s, t) == a*DiracDelta(t - 3)
    assert ILT(exp(2*s), s, t) == InverseLaplaceTransform(exp(2*s), s, t, None)
    r = Symbol('r', real=True)
    assert ILT(exp(r*s), s, t) == InverseLaplaceTransform(exp(r*s), s, t, None)
    # Rules for testing whether Heaviside(t) is treated properly in diff rule
    assert ILT(s**2/(a**2 + s**2), s, t) == (
        -a*sin(a*t)*Heaviside(t) + DiracDelta(t))
    assert ILT(s**2*(f(s) + 1/(a**2 + s**2)), s, t) == (
        -a*sin(a*t)*Heaviside(t) + DiracDelta(t) +
        Derivative(InverseLaplaceTransform(f(s), s, t, None), (t, 2)))
    # Rules from the previous test_inverse_laplace_transform_delta_cond():
    assert (ILT(exp(r*s), s, t, noconds=False) ==
            (InverseLaplaceTransform(exp(r*s), s, t, None), True))
    # inversion does not exist: verify it doesn't evaluate to DiracDelta
    for z in (Symbol('z', extended_real=False),
              Symbol('z', imaginary=True, zero=False)):
        f = ILT(exp(z*s), s, t, noconds=False)
        f = f[0] if isinstance(f, tuple) else f
        assert f.func != DiracDelta


@slow
def test_expint():
    x = Symbol('x')
    a = Symbol('a')
    u = Symbol('u', polar=True)

    # TODO LT of Si, Shi, Chi is a mess ...
    assert laplace_transform(Ci(x), x, s) == (-log(1 + s**2)/2/s, 0, True)
    assert (laplace_transform(expint(a, x), x, s, simplify=True) ==
            (lerchphi(s*exp_polar(I*pi), 1, a), 0, re(a) > S.Zero))
    assert (laplace_transform(expint(1, x), x, s, simplify=True) ==
            (log(s + 1)/s, 0, True))
    assert (laplace_transform(expint(2, x), x, s, simplify=True) ==
            ((s - log(s + 1))/s**2, 0, True))
    assert (inverse_laplace_transform(-log(1 + s**2)/2/s, s, u).expand() ==
            Heaviside(u)*Ci(u))
    assert (
        inverse_laplace_transform(log(s + 1)/s, s, x,
                                  simplify=True).rewrite(expint) ==
        Heaviside(x)*E1(x))
    assert (
        inverse_laplace_transform(
            (s - log(s + 1))/s**2, s, x,
            simplify=True).rewrite(expint).expand() ==
        (expint(2, x)*Heaviside(x)).rewrite(Ei).rewrite(expint).expand())
